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We investigate, in three spatial dimensions, the transition from the normal state to the Fulde- 
Ferrel-Larkin-Ovchinnikov superfiuid phases. We make use of a Fourier expansion for the order 
parameter and the Green's functions to handle the quasiclassical equations in the vicinity of the 
transition. We show that, below the tricritical point, the transition is always first order. We find 
that, at the transition, the higher Fourier components in the order parameter are always essentially 
negligible. Below the tricritical point we have the already known result that the order parameter 
has a spatial dependence which is essentially cos(q.r). However when the temperature is lowered, 
the order parameter switches to a sum of two cosines, with equal weigths and wavevector with the 
same length, but orthogonal directions. Finally by further lowering the temperature, and down to 
T = 0, one finds a another transition toward an order parameter which is the sum of three cosines 
with again equal weigths and orthogonal directions. Hence the structure of the order parameter gets 
more complex as the temperature is lowered. On the other hand the resulting critical temperatures 
are found to be only slightly higher than the ones corresponding to the standard second order 
FFLO transition. We apply our results to the specific case of ultracold Fermi gases and show that 
the differences in atomic populations of the two hyperfine states involved in the BCS condensation 
display sizeable variations when one goes from the normal state to the superfiuid FFLO phases, or 
one FFLO phase to another. Experimentally this should allow to identify clearly the various phase 
transitions. 

PACS numbers : 05.30.Fk, 67.90. +z, 74.20.Fg, 74.25. Op 



I. INTRODUCTION 

Although initiated a long time ago [1,2] the problem of the Fulde-Ferrel-Larkin-Ovchinnikov phases is far from 
having found a complete solution. In its simplest version one disregards the standard coupling of the magnetic 
field to the orbital degrees of freedom and one looks at the coupling with the electronic spins which gives rise to a 
difference between the chemical potentials of the up and the down spin populations. The question is to know how the 
superconducting transition is modified under these circumstances, and naturally the answer is relevant to deal with the 
more complex structures expected to appear when the orbital coupling is taken into account and vortex-like structures 
will be produced [3] . A recent and quite interesting development has been the realization that this FFLO phases are 
equally of high interest for the physics of ultracold atomic gases in their superfiuid BCS-like state [4] , and also in the 
physics of quark matter and dense nuclear matter [5-7] as it might be found in the center of neutron stars. Our results 
near the transition at higher temperature [8] have confirmed the earlier results found in the literature [9-11]. Namely 
the preferred order parameter is a one-dimensional order parameter, whether the space itself is three-dimensional (3D) 
or two-dimensional (2D) (as it is relevant for lamellar superconductors or cuprate superconductors). At the transition 
it is proportional to a simple cosine A(r) ~ cos(q.r), this result being exact in 2D and almost exact in 3D. However 
we have found recently [12] that, in 2D, the situation gets increasingly complex at low temperature. Indeed in this 
case the order parameter at the transition takes the form of a superposition of cosines, with wavevectors having the 
same length but different directions, the number of these cosines increasing indefinitely when the temperature is going 
to zero. 

In the present paper our purpose is to explore the same problem in three dimensions. This question is technically 
much more difficult than in 2D, since the transition turns out to be always first order in contrast with the 2D case. 
This fact is already known from higher temperature results [9-11,8] and we will find that it remains true down to 
T = 0. This makes it impossible to make use of a Ginzburg-Landau type of expansion, which in 3D is valid only 
near the tricritical point and, in 2D, at any temperature. One is thus forced, even at the transition, to face the 
full non-linear problem, which is quite difficult to handle. A most convenient way to attack it is to make use of the 
quasiclassical equations of Eilenberger [13], Larkin and Ovchinnikov [14]. Nevertheless these equations themselves are 
not so easy to manipulate, even numerically, with the prospect of finding possible solutions with a three-dimensional 
order parameter, in a way analogous to what we have done in 2D. In a preceding paper [15] we have made a progress in 
this direction by showing that the introduction of a Fourier expansion in the quasiclassical equations allows to obtain a 
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solution which converges very rapidly toward the exact result. As a consequence a few terms in the expansion provide 
an excellent approximation. In this preceding paper only the principle of this method and its application to the a 
transition with one-dimensional order parameter have been considered. It has been checked against results obtained 
by Matsuo et al [10], which started also from the quasiclassical equations, but used another method to solve them. We 
have found excellent agreement. In the present paper we will apply our technique to more complex cases where the 
order parameter does not have anymore a simple one-dimensional dependence. Naturally this was actually the aim 
of introducing this Fourier expansion method. We will explore which order parameter is actually favored by FFLO 
phases, in three dimensions, at low temperature. We find that this is no longer a one-dimensional order parameter, 
as it found at higher temperature, with a spatial which is nearly given by A(r) ~ cos(q.r) at the transition. We will 
rather obtain an order parameter with a more complicated structure. At the transition it is nearly the sum of two or 
three cosines (this last case being found near T = 0), with directions of oscillations orthogonal to each other. As a 
result the transition does no longer switch to a second order one at low temperature, as found by Matsuo et al [10], 
and go to the second order transition originally investigated at T = by Larkin and Ovchinnikov [2] . It stays rather 
a first order transition down to T = 0. Our results have been briefly reported elsewhere [16]. 

Finally we apply our results to the case of ultracold Fermi gases and calculate the difference between the atomic 
populations of the two hyperfine states involved in the formation of Cooper pairs. We find that these differences can 
be quite sizeable and allow experimentally a clear distinction between the various phases, normal or superfluid, which 
may be found in these systems. In particular this should allow experimentally a clear identification of the transitions 
between these phases. 

II. FORMALISM 

In this section we will briefly recall for completeness the main points of the method presented in Ref. [15], hereafter 
referred to as I. However we will write it immediately for the most general case we have in mind, i.e. a general Fourier 
expansion for the order parameter A(r) and the Green's functions. More explicit formulations can be found in I. 

We write Eilenberger's equations for the diagonal g(u),\t, r) and off-diagonal /(u>,k, r) quasiclassical propagators. 
Here r is the spatial variable while k is the angular location of the wavevector on the spherical Fermi surface of 
radius hp. Finally lo is a general Matsubara frequency, to be continued into u> n — ip,, where u> n = (2n + l)wT is the 
standard Matsubara frequency and 2/2 = /ij — /ij is the chemical potential difference between spin up and spin down 
populations of the particles forming Cooper pairs. These equations read [13]: 

(w + k.V)/(w,k,r) = A(r)</(w,k,r) 
(w - k. V)/ + (w, k, r) = A* (r)g(u, k, r) 
2k.Vff(w,k,r) = A*(r)/(a;,k,r)-A(r)/ + (a;,k,r) (1) 

in which we have simplified the writing by taking h — 1 and m = 1/2 for the particle mass. These linear equations 
are closed by the normalization condition : 

g{u>, k, r) = (1 - f(oj, k, r)/+ ( w , k, r)) 1 ^ (2) 

from which one can see that the last of Eq.(l) results from the two first ones. Finally the order parameter itself is 
linked to /(w,k, r) by the self-consistency equation [13,15] which we will not write explicitely here. 
We assume then a very general Fourier expansion form for the order parameter: 

A(r) = A {«a ex P(* m<\i*) (3) 

{«.} 

where the rii are relative integers with i — 1,2, ...,N q and the qi are N q general wavevectors for which no specific 
relations are assumed. If we had only three wavevectors with i = 1,2,3 this would mean that we consider a general 
periodic order parameter. However we do not restrict ourselves to three wavevectors, and we will indeed consider 
below explicitely the case of four wavevectors. We note that such an order parameter is not periodic in general. 
Naturally the numerics gets rapidly more difficult when the number of wavevectors qi increases. We will nevertheless 
assume for simplicity a real order parameter since, at least in three spatial dimensions, all the order parameters found 
for actual FFLO phases have this property. This implies At, = A{_ ni j. We also restrict ourselves as before to an 
order parameter having real components A{„.}. Finally we will only explore the case where all the wavevectors have 
equal length |qj| = q. Indeed we know from the original work [2] of Larkin and Ovchinnikov that the physics of the 
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FFLO phases is that there is a competition between the different directions of the wavevectors having same modulus. 
At least this is the prevalent situation near the transition, in which we are interested in. 
We make for the Green's functions similar expansions: 

/( r ) = X) •/W ex p(*X niCti - r ) / + ( r ) = E{„ i }/{k i } ex p( i S i «iqi- r ) g( r ) = X^{"i} ex p(*X] niqi - r ) ( 4 ) 

{m} i {rii} i 

and substitute them into Eilenberger's equations Eq.(l). We introduce for convenience d{ ni } = (/{«;} — 
which gives f{ ni } = {i — W / K ({ n i}))^{n;}j where n({ni}) — k. J^i n i c Li- This leads to: 

n{{ni\) ^ 

{ "* } = ~ LU 2 + K 2 ({n,|) ^ A {Pi}W{ni- Pi } + .9{r ll+Pl }) 

■ 9{ ™ l} = Ac({nJ) ^ A {p.}( d {™»-P.} + ( 5 ) 

Just as in I, we consider only order parameters with zero spatial average i.e. A{ } = when all the rii = 0, because 
mixing in the uniform BCS phase order parameter is likely to be unfavorable. Then it can be seen iteratively, as in 
I, that only odd components of d^ n .j and even components of g{ ni } arc non zero, i.e. precisely = if ^ i n, is 

even, and g{ ni } = if ^ rii is odd. 

Eq.(5) look quite complicated, but it turned out that we had to deal with them only in rather simple cases. Indeed 
we studied first the case of two cosines, and concluded as in I that the weight of higher order harmonics for the order 
parameter is quite small (we will give specific results below). We have then taken this conclusion for granted when we 
studied three or four cosines. This means that we have mostly dealed with Eq.(5) only in the case where all the A{„. -j. 
are zero, except when ^\ |rjj| = 1, i.e. all the rii are zero except one of them which is ±1. In this case we have only 
2N q terms in the sums in the right-hand side of Eq.(5). Similarly we have studied in some cases for two cosines the 
situation where the weights A{ ni } are not equal and we have come with the conclusion that the symmetric situation 
with equal weights is the most favorable. We have then taken this conclusion for granted in the rest of our study. This 
means finally that, in all these cases, we have to deal with a single non zero component of A/„.\, which we call Ai in 
the following. Naturally, apart from the complexity of the equations, the essential difficulty in this FFLO problem is 
that the functional space to be explored is huge, to say the least, and that it seems hopeless to explore it completely. 
One has thus to appeal to symmetry arguments, like the equal weight argument, to restrict it. 

Our two sets of Fourier components, g{ ni y and c?{ ni }, are defined on a discrete A^-dimensional lattice. However 
since on any site one of these two components is zero depending on the parity of J2i n i> as we have seen, we have 
actually to deal on each site with a single component. Then Eq.(5) relate the value of this component on a given site 
to its a value on the nearest neighbours of this site on the 7V g -dimensional lattice. By the same kind of arguments 
as in I it can be seen that the Fourier components decrease very rapidly on the site goes far away from the origin. 
Hence we will find an approximate solution for our problem by requiring that our components are zero when the site 
is far away from the origin. Specifically we will take them to be zero when \rii\ > N max . This serves as boundary 
condition for our linear system Eq. (5) . In this way equations become a finite dimensional linear system, which we solve 
numerically by standard and efficient methods. In principle when we let N max go to oo we find the exact solution. 
In practice, as will be discussed below, the convergence is found to fast enough so that the size of our linear system 
stays quite reasonable. We note that, because we have the properties g{- ni } = <?{«;} an d ^{-ni} = ~d{ ni y, we have 
only to deal with half of the sites satisfying ^ \m\ < N max . This domain T> for summation over the {rii} is shown 
cxplicitcly in Fig. 1 in the case of two cosines N q = 2 and for N max = 3 . 
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FIG. 1. Lattice sites, in the {ni,ri2} space, involved in the summation in Eq.(7), in the case of two cosines N q = 2 and 
for N m ax = 3 . The numbers correspond to the actual ordering we use in writing the corresponding equations, which goes by 
increasing |ni| + |ri2 [- 

In practice we have set, as in I: 

d{ n ,} = K{{ni})D {nz} g 

9M = Ginkgo (6) 
where go is for g{ ni= o}, i- e - the value of g{ ni } at the origin of our discrete lattice. Then Eq.(5) becomes explicitely: 

[u 2 + K 2 ({m})}D {m} + A x G M = 
- K ({ ni })G {ni} + <{nj})D {nj} = (7) 

{ nj } 

We normalize these new components by the condition Go = 1. When the normalization condition Eq.(2) is expressed 
with these components, one finds: 

. 9fi 7 2 = 1 + 2^Gf ni} + 2^> 2 - K 2 {{ ni })]D\ ni} (8) 
v v 

This component g , the spatial average of the Green's function, is just the quantity we need to calculate [15] the free 
energy difference Q s — VL n between the superfluid and the normal phase: 



q _ o t °° r°° r m a 2 



where No is the single spin density of states at the Fermi surface. In Eq.(9) the angular average is over the direction 
of k at the Fermi surface Sf- We have also, instead of the standard BCS critical temperature T c o at zero effective field 
jl = 0, introduced the critical temperature T sp (p,/T) for the spinodal transition, which is the second order transition 
from the normal to the uniform BCS phase. This T sp is conveniently seen as a function of p,/T, which is a kind of 
angular coordinate when one explores the (jl, T) plane. More specifically, since we are only interested in the location 
T c of the first order transition, we look at fixed fi/T for the highest temperature T = T c at which we will have 
Q s — Q n = 0. Naturally we maximize this temperature with respect to the amplitude Ai of the order parameter and 
to the common length q of the wavevector q,; . 

III. NUMERICAL RESULTS 

We start the presentation of our numerical results by the case of two wavevectors N q = 2. Indeed as already 
mentionned we have in this case made specific explorations, the results of which we have taken for granted in the 
numerically more heavy situations of three or four wavevectors. 
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A. Two wavevectors 



We have investigated how the critical temperature T depends on the angle between the two wavevectors qi and q 2 . 
In this study we have assumed that only the lowest Fourier components arc important, so that our order parameter 
is actually the sum of two cosines corresponding to these wavevectors. Our results are given in Fig. 2. As in I, we 
plot the ratio of the critical temperature T/Tfflo to the FFLO critical temperature Tfflo{P>/T) obtained for the 
same value of the ratio p/T. Instead of fi/T, we give on the x-axis the value of TpFLoifi/T) itself, compared to the 
standard BCS critical temperature T c q. 
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FIG. 2. Critical temperature T for the order parameter A(r) = 2Ai[cos(qi • r) + cos(q2 • r)], for a selection of angles a 
between the wavevectors qi and q2- Precisely we plot the ratio of the critical temperature T/Tfflo to the FFLO critical 
temperature Tfflo{P-/T) obtained for the same value of the ratio jX/T. On the x-axis, instead of fl/T, we give the value of 
Tfflo(P,/T) itself, compared to the standard BCS critical temperature T c o- 
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We give results only for selected angles between and 7r/2, since the results are the same for a and 7r — a. We also 
give the result for a = which is just the result for a single cosine. It can be seen that, above Tfflo/T c o ~ 0.154, 
this is this single cosine order parameter which is favored. Below this temperature the order parameter switches 
discontinuously to two cosines, with an angle a = tt/2 between the corresponding wavevectors. Hence there is no 
situation where an angle < a < n/2 is favored. We note that the limit a — > is singular. This can be seen on 
Fig. 2 where this limit is noted a/0. This singularity can be understood since, for two cosines, the spatial average 
of the square of the order parameter is always 4A^, even if the angle between the two vectors is very small, while it 
is 8A^ when this angle is exactly zero. In the following we take this symmetry result for granted and assume that 
the situation with orthogonal wavevectors is always the most favorable. We have also made some exploration of the 
evolution of the critical temperature when the relative weights Aj of the cosines vary. For example, we have kept the 
wavevectors orthogonal or with an angle n/3 and minimized independently the two weights. We have found that, in 
either case, the optimum is for equal weights. 

We consider next the sensitivity to our upper bound N max to the order of Fourier components we take into account, 
that is we study how fast our procedure converges to the exact result, obtained in principle for N max = 00. We make 
this study only for the most favorable situation a — n/2, which is the only one we consider from now on. The result 
is displayed in Fig. 3, where we plot, as a function of the temperature, the results for N max = 3,4 and 5. The results 
for N max — 4 and N max — 5 are undistinguishablc within our precision, and it can be seen that N max — 3 gives 
already an excellent precision. This is actually better than what we have found in I, for the case of a single cosine. On 
the other hand the numerics is more complicated in the case of two cosines because the angular average requires two 
angular integrations, one on the polar angle and the other one on the azimuthal angle, while for symmetry reasons 
only the polar integration is necessary for a single cosine. We note that, in this respect, the situation does not get 
worse when we consider three or four cosines, since we have still only two angular integrations to perform. 
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FIG. 3. Critical temperature T for the order parameter A(r) = 2Ai[cos(qi ■ r) + cos(q2 ■ r)], for an angle a = tt/2 between 
the wavevectors qi and q2. The full line is for N max = 3, while the filled circles and the open diamonds are for N max = 4 and 
5. Actually these last results Nmax = 4 and 5 can not be distinguished within our precision. 




We also give for completeness in Fig. 4 the corresponding results for the size Ai//2 of the order parameter and the 
corresponding length q = qkF/(2rnp) of the wavevectors. The results are somewhat scattered, since our minimization 
procedure makes it difficult to find exactly the location of the minimum. On the other hand the critical temperature 
is much more precisely known. 
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FIG. 4. Relative amplitude Ai/p. of the order parameter and corresponding length q = qltF / (2mjl) of the wavevectors, as 
a function of the upper bound N max for the order of Fourier components. The full line is for N max = 3, the filled circles for 
Nmax = 4 and the open diamonds for N ma x ~ 5. 
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Finally we consider how our results are modified when we go beyond the lowest harmonic approximation. We take 
for this purpose the order parameter: 



A(r) = 2Ai[cos(qi • r) + cos(q 2 • r)] + 2A 3 [cos(3qi • r) + cos(3q 2 • r) 



(10) 



and study the size of the higher harmonic correction A3 as a function of temperature. This is done for Nmax=4. The 
result is found in Fig. 5. We see that the harmonic A3 is always quite small, since the ratio A3/A1 is typically a few 
10 -3 . Hence, just as in the case of the one-dimensional order parameter [15] , it is quite justified to consider only the 
lowest harmonic. 




FIG. 5. Relative amplitude of the higher harmonic A3//X for the order parameter Eq.10. 

We give also in Fig. 6 the results for the critical temperature, the lowest harmonic amplitude and the wavevector 
length with and without taking into account the higher harmonic A3. One can see that the results are indeed 
essentially undistinguishable. 
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FIG. 6. Critical temperature, lowest harmonic amplitude Ai / fi and corresponding wavevector length q — qk F / '(2mfi) for the 
order parameter Eq.10 with A3 = and A3 / (as given by Fig. 5) respectively. 



B. Three wavevectors 



We consider now the case of three wavevectors qi, q2 and q3. Following the above study, we take for granted that 
higher harmonics give a negligible contribution so that we can restrict ourselves to the order parameter: 

A(r) = 2Ai [cos(qi • r) + cos(q 2 • r) + cos(q 3 • r)] (11) 

which implies that we also assume the amplitude for the three cosines are equal. Similarly we assume that the optimum 
order parameter corresponds to wavevectors with same length |qi| = | q2 1 = | C3.3 1 - Finally we have again studied in 
this case how the critical temperature depends on the angle between the three wavevectors, by taking them in a 
rhomboedral geometry and varying the rhomboedric angle. We have found again that the most favorable situation is 
found when the wavevectors are orthogonal. In the same way as for two wavevectors, this result is consistent with a 
phcnomcnological interpretation in terms of an effective repulsive potential between any two wavevectors, decreasing 
when the angle between them increases. With such a description the total potential is clearly minimized when the 
three wavevectors are orthogonal. 

We present now our numerical results for this specific situation of orthogonal wavevectors. The critical temperature 
is given in Fig. 7. We display again for comparison on this figure the results for a single cosine and for two cosines, 
given in Fig. 2. We see from this figure that the order parameter at the transition switches from a single cosine 
for T FFL o/T c0 > 0.154 to an order parameter with two cosines for 0.080 < T FFLO /T c0 < 0.154, and to an order 
parameter with three cosines for T FFF o/T c o < 0.080. Naturally, since they do not have the same symmetry, this 
implies that, inside the superfluid phase, there are first order transition lines between these various order parameters 
(or rather their continuation, because higher harmonics are expected to be more important deeper in the superfluid 
phase). On the other hand it is rather striking that the switch of the order parameter from a given symmetry to 



8 



another has very little effect on the value of the critical temperature itself since T c /Tfflo is only slightly larger than 
unity. 
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FIG. 7. Critical temperature for an order parameter with respectively three cosines (short dashed), two cosines (long dashed) 
and a single cosine. The actual transition between the normal and the superfluid state corresponds to the highest possible 
transition temperature between these three possibilities. The filled dots refer to the "cube" configuration, considered in the 
next subsection. 




The results in Fig. 7 for three cosines have been obtained for N max = 3. We have made several checks, by performing 
calculations for N max = 4, that convergence is indeed already reached. This is shown in Fig. 8 where the results for 
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FIG. 8. Critical temperature for an order parameter with three cosines. The full line is for N max 
are for N max = 4. 
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Finally we present for completeness in Fig. 9 our corresponding results for the relative amplitude Ai//z and the 
reduced wavevector length q of the order parameter. 
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FIG. 9. Order parameter amplitude Ai//i and corresponding reduced wavevector length q as a function of temperature. 

C. Four wavevectors 

As we have seen above, when T — > 0, increasing the number of cosines from a single cosine to three increases the 
critical temperature. It is natural to wonder if this trend does not keep going for a larger number of cosines, all the 
more since this is indeed what happens in two dimensions [12] when T — > 0. In order to explore this possibility we 
have considered the transition toward an order parameter which is the sum of four cosines: 

A(r) = 2Ai [cos(qi • r) + cos(q 2 • r) + cos(q 3 • r) + cos(q 4 • r)] (12) 

with same wavevector length |qi| = | q2 1 = |q.3| = |q4|- We have considered only the most symmetri- 
cal situation, where these wavevectors point toward the corners of a cube, corresponding to the directions 
±(1, — 1, — 1), ±(— 1, 1, — 1), ±(— 1, — 1, 1) and ±(1,1,1). The angle between any of these wavevectors is 70.5°. We 
have restricted our calculations to the case N max — 3. Our results are given in Fig. 7, as the "cube" configuration. It 
is clear that this configuration is never favored. Hence when the temperature goes to zero the optimal configuration 
has three cosines, with wavevectors pointing in orthogonal directions. We discuss below a possible phenomenological 
interpretation of this result. 

D. Approximate N max — 1 solutions 

Since it is not so easy to find an insight into the complex solutions obtained numerically, it is worthwhile to consider 
the simplest meaningful approximation obtained by retaining only the lowest order Fourier components of the Green's 
function, that is taking N max = 1. A major advantage is that it can be handled easily analytically, which could be 
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useful to reach a deeper understanding of the solutions for the FFLO phases. For example the n = component of 
the diagonal propagator is given by: 



9o 



(k ■ q; 



-1/2 



[CJ 2 + (k • q; 



(13) 



for a general order parameter, which is the sum of n cosines with weights Aj and wavevectors q.;. Naturally this 
makes also the final numerical treatment much easier. 

In the case of a single cosine we have seen [15] that this solution gives the proper qualitative behaviour, with a 
switch from a first order to a second order transition, although the location of the switching temperature is not very 
accurately given. We perform here the same kind of study in the case of two, three and four cosines. In Fig. 10 we 
compare the results of this N max = 1 solution with our complete solution, presented in the preceding subsections. 



A(r) = 2Ao \cos(qx) + cosiqy) + cos(qz) 
A(r) = 2Aq [cos(qx) + cos(qy) 
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FIG. 10. Upper panel: critical temperature for our two cosines order parameter with N ma x = 1 approximation (dotted line) 
compared to our complete result (dashed-dotted line). Same comparison for our three cosines order parameter : dashed line 
(N m ax = 1 approximation), full line (complete result). Lower panel : same comparison for our four cosines order parameter : 
dotted line (N max = 1 approximation), filled circles (complete result). 



It can be seen that this N max — 1 approximation is quite satisfactory since it reproduces quite well qualitatively 
and semiquantitativcly our complete results. Nevertheless this is compounded with the fact that all the results are 
quite near unity. Hence small differences, which could by themselves be considered as unsignificant, can lead to 
qualitative differences in the final results. This is displayed in Fig. 11 where we gather the results for the various 
order parameters within the N max — 1 approximation. It can be seen that the three cosines order parameter always 
dominates at low temperature, in contrast to our complete results where there is a temperature range where the two 
cosines order parameter is the best one. Similarly the four cosines order parameter is better than the two cosines one 
at very low temperature within the N max = 1 approximation while this never occurs with our complete solution. In 
conclusion this N max — 1 approximation is quite convenient for a fast exploration of the free energy minimization 



11 



problem. It gives a semiquantitatively correct picture of the competition between various order parameters. However 
it can not be fully trusted quantitatively. 




E. Phenomenological interpretation 

We compare now our results to the phenomenological picture suggested by Bowers and Rajagopal (BR) [5]. In their 
T = study BR made use of a Ginzburg-Landau expansion to find the most favored order parameters in the case of 
a second order transition, even if in many cases the transition is actually first order which invalidates the expansion. 
They came in the course of this work to the following heuristic picture. To each wavevector present in the order 
parameter, one associates a circle drawn on the unit sphere, the axis of this circle being the wavevector direction and 
its angular opening ip being given by cos("0o/2) = 1/? where q ~ 1.2 is the reduced wavevector [1,2] for the T = 
second order FFLO phase transition. This gives ipo ~ 67.1°. This circle, in the Fulde-Ferrell description, corresponds 
to a favored region for pairing (it is infinitely thin because we are just at the second order phase transition where the 
order parameter is zero). A first principle in BR heuristic picture is that the most favored order parameter has the 
maximum number of wavevectors, in order to increase the number of circles, i.e. the domains where favorable pairing 
occurs. On the other the second principle is that the crossing of two circles is very unfavorable energetically and must 
be avoided. This no-crossing principle of circles with a definite angular size leads clearly to an upper bound in the 
number of possible wavevectors. It has been found by BR to be 9 wavevectors. However the corresponding situation is 
not a symmetrical one for the wavevectors, and accordingly BR suggested that the symmetrical form, corresponding 
actually to our 'cube' configuration, would be the most favorable. This can be understood if the no-crossing principle 
is considered as a manifestation of an effective repulsive interaction between wavevectors directions. In this case it is 
reasonable to believe that the repulsion will be minimized for a symmetrical geometry for the wavevectors. 

Now we have just found that the most favorable order parameter at T = has 6 wavevectors (corresponding to 
our 3 cosines order parameter) instead of 8. But this can be understood qualitatively as resulting from the fact 
that the transition is actually first order, instead of second order, which makes the order parameter non zero at the 
transition. In such a case, instead of being infinitely thin, the region corresponding to favorable pairing will widen 
with a thickness linked to the non zero value of the order parameter. This will make the circles effectively larger, 
which may explain why it is no longer possible to fit 8 circles on the unit sphere and it is necessary to reduce their 
number to 6. Hence there is a reasonable agreement between this phenomenology and our results. We note finally 
that the results we have found in two dimensions [12] support also this picture. Here the circles on a unit sphere are 
replaced by to points on a unit circle. In this case the transition is second order, and as T — > the distance between 
the two points goes to zero. This allows to fit on the unit circle an ever increasing number of wavevectors, as T — > 0. 
This is indeed what we have found [12] . 



IV. THE CASE OF ULTRACOLD FERMI GASES 



We will now apply our results to the case of ultracold gases, which is presently a field of very strong interest both 
experimentally and theoretically. So far our analysis has been restricted to the case where the chemical potential 
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difference between the two spin species is fixed. This is of course the reievant situation for superconductors where 
this chemical potential difference is produced by very high magnetic fields, and the resulting phase diagram which we 
have found above applies directly to this case. On the other hand in ultracold gases these are rather the populations 
of the various hypcrfine states of the atom under study which are under direct experimental control. Accordingly, in 
order to apply our results to the case of ultracold atomic Fermi gases, we need to investigate more specifically the 
effect of a fixed density difference on the phase diagram, rather than a chemical potential difference. 

This modification is quite relevant in our 3D case since we have essentially to deal with first order transitions. This is 
in contrast with the 2D case where the transitions are found [12] to be second order. Across a second order transition, 
physical quantities like the order parameter or the density are continuous. This implies that phase diagrams for 
fixed chemical potentials or fixed densities are essentially equivalent. On the other hand the order parameter and the 
density are discontinuous for a first order transition, as for example in the well known case of the standard liquid-gas 
transition. There is a forbidden domain in the phase diagram, corresponding to phase separation, where the two 
phases, normal and superfluid, are coexisting in the system. This amounts, so to speak, to split the transition line 
which appears for fixed chemical potential difference into two lines in the case of fixed density difference. Each of 
this line gives the density of one of the coexisting phases. In our case the situation is somewhat more complicated 
since depending on the temperature we have different phases coming into play. Hence we consider now the density 
difference in the various FFLO phases we have found. 

From the solutions of the Eilenberger equations that we have obtained using our Fourier expansion method, it is 
quite straightforward to calculate any physical quantities. In particular the difference in the densities for spin up and 
down is given by 



"t( r ) ~"i( r ) _ 9l - 
N - 2fi 



47rTIm V" < g(ui n - ifj,, k, r) - g n (uj n - ifj, k) > jc 



=o 



(14) 



where the first term in the right-hand side, 2j2, is actually the result for the normal state and g n is the normal state 
Green's function. It is worth noting that this density difference shows oscillations at spatial frequencies equal to 2q, 
Aq, and so on, which we have found numerically to be quite sizeable. These oscillations could be used as signatures of 
the FFLO phases, in analogy with the oscillations of the magnetization in the case of superconductors. However we 
will restrict ourselves to macroscopic physical quantities, and accordingly we consider only the average value of the 
density difference which amounts to replace g in Eq. (14) by go. 

The calculation of this average density difference for the different relevant phases leads us to the phase diagram 
shown in Fig. 12. 



T/T c ° 0.5 




0.5 0.55 0.6 0.65 0.7 0.75 0.8 

An/(2N A ) 

FIG. 12. Phase diagram for reduced temperature T/T c o versus atomic population difference An. The transition lines for 
cases 1, 2 and 3 (see text) are respectively represented by full lines, by full circles and by diamonds. The dashed-dotted line 
represents the standard second order FFLO transition which goes into the standard BCS transition above the tricritical point. 



Just as in the previous case where the chemical potential was fixed, we have three different temperature domains 
below the tricritical point: 

1. The region 0.154 T c0 < T < T tcp ~ 0.561 T c0 where the transition is from the normal state to the superfluid 
state with A(r) = Ai cos{qx). 
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2. The region 0.08 T c0 < T < 0.154 T c0 where one goes from the normal state to the A(r) = Ai (cos(qx) +cos(qy)) 
state. 

3. The region T < 0.08 T M where the order parameter of the superfluid state is A(r) = Ai(cos(gx) + cos(qy) + 
cos(qx)). 

The most remarkable feature of these results is that the differences for the value of An between the normal and the 
various FFLO phases are quite sizeable, which should allow experimentally to locate clearly the BCS transition. The 
other important point is that, between the various FFLO phases themselves, An changes in a quite noticeable way. 
Its evolution is qualitatively easy to understand. The important differences in atomic populations come [9] from the 
regions where the order parameter is small while, in the regions where it is large, pairing favors more equal populations. 
Since the number of cosines is increased when we go to lower temperatures, this implies that the regions with small 
order parameter are more important. This is coherent with the fact that, at low temperature, the superfluid phase 
allows a larger difference between atomic populations than the superfluid phases found at higher temperature. This 
feature should allow to locate rather easily the change of structure of the order parameter, while as we have seen 
the changes in the critical temperature between the various phases is quite small and would not allow an convenient 
observation. 



V. CONCLUSION 



In this paper we have considered the original problem raised by Fuldc, Ferrell, Larkin and Ovchinnikov, namely the 
possibility and the nature of the transition from a normal Fermi liquid to a superfluid phase with a space dependent 
order parameter, when the chemical potential of the two spin populations are different. The Fermi surface is assumed 
to be spherical, in three spatial dimensions. As the transition turns out to be first order we have made use of the 
quasiclassical Eilenberger's equations to handle this non linear problem. In order to solve these equations we have 
introduced a Fourier expansion for the order parameter as well as for the Green's functions. In a preceding paper [15] 
we have introduced in details this procedure, we have checked that it worked properly on the case of a one-dimensional 
order parameter and that it is quite efficient, at least in the vicinity of the transition. 

In the present paper we have made use of it to show that, at low temperature the order parameter switches from the 
one-dimensional form, found at higher temperature, which at the transition, is essentially a simple cosine, to a more 
complex structure. At zero temperature we find that the most stable order parameter is essentially, at the transition, 
the sum of three cosines with equal weights and orthogonal wavevectors. When the temperature is raised, the order 
parameter switches to a sum of two cosines, also with equal weights and orthogonal wavevectors, before reducing to a 
single cosine at higher temperature. Hence we have obtained that, at the transition, the higher Fourier components in 
the order parameter are essentially negligible. This makes clearly our procedure particularly efficient. We have found 
that, actually, the critical temperatures for the first order transitions toward these phases is only slightly higher than 
the critical temperature for the standard second order FFLO transition. This seems to show that there are many 
phases with nearly equal free energy. This should be for example included in the theory when, for superconductors, 
the coupling of the magnetic field to orbital degrees of freedom is also taken into account. We stress again that our 
method provides a way to obtain some insight in the analytical structure of this complex theory. Finally we have 
applied our results to the specific case of ultracold fermionic atoms and shown that the various FFLO phases arising 
at low temperature give rise to differences in atomic populations which allows to identify them easily. 

* Laboratoire associe au Centre National de la Recherche Scientifique et aux Universites Paris 6 et Paris 7. 
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